package lib;

import java.util.HashMap;

/**Neighborhood Correlation common routines
 * (see http://www.neighborhoodcorrelation.org)
 *
 * (C) 2011 Jacob Joseph <jmjoseph@andrew.cmu.edu>
 *          and Carnegie Mellon University
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation; either version 3
 * of the License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 */
public class NC_base {
    HashMap<String, String> blast_info;          // dict of BLAST run parameters
    HashMap<String, Double> nc_info;             // dict of NC run parameters    
    
    HashMap<String, String> var_cache;           // Sequence variance: {'seq_id': var_x, ...}
    HashMap<String, Triple<String, String, String>> E_cache;             // {'seq_id': (n_x, sum_x, E_x), ...}

    public NC_base() {
    	blast_info = new HashMap<String, String>();
        nc_info = new HashMap<String, Double>();

        var_cache = new HashMap<String, String>();
        E_cache = new HashMap<String, Triple<String, String, String>>();
    	
 /*   	NChelparr nchelparr = new nchelparr();
    	nchelparr.calculate_nc_inner();
  */  }
    
    public void new_run(String br_id, String e_thresh, String nc_thresh,
    		String smin_factor, String self_hits, String num_residues,
    		String score_type) {

    	HashMap<String, String> nc_info = new HashMap<String, String>();
    	nc_info.put("br_id", br_id);
    	nc_info.put("e_thresh", e_thresh);
    	nc_info.put("nc_thresh", nc_thresh);
    	nc_info.put("smin_factor", smin_factor);
    	nc_info.put("score_type", score_type);
    	nc_info.put("self_hits", self_hits);

    	if (num_residues != null)
    		blast_info.put("num_residues", num_residues);
    }

    public void calculate_smin() {

    	assert blast_info.size() > 0; //blast_info not initialized. fetch_score_dict_flatfile() must first be executed, or a database call made.

    	// Fetch blast run information
/*    	Double num_residues = Double.parseDouble(blast_info.get("num_residues"));
    	Double blast_e_param = blast_info.get("params"]["expectation"];
   		Double smin_factor = nc_info.get("smin_factor");

   		double smin = (Math.log( (Math.pow(num_residues, 2)) / blast_e_param) / Math.log(2));

   		// reduce the true smin by some factor
   		smin = smin * smin_factor;
   		nc_info.put("smin", smin);
   		nc_info.put("logsmin", Math.log10(smin));*/
    }

/*    public Triple<String, String, String> E_x(HashMap<String, Couple<String, Double>> score_dict, String seq_id) {
        //Calculate and cache n_x, sum_x, and E_x

        if (E_cache.containsKey(seq_id))
            return E_cache.get(seq_id);

        if(!nc_info.containsKey("smin"))
        	calculate_smin();
        Double scores = score_dict.get(seq_id).second;

        int num_seqs = Integer.parseInt(blast_info.get("num_sequences"));
        logsmin = nc_info.get("logsmin");

        n_x = scores.size;
        others = num_seqs - n_x;
        Double sum_x = Math.sum(scores);
        E_x = (sum_x + logsmin * others) / num_seqs;

        E_cache[seq_id] = (n_x, sum_x, E_x);
        return (n_x, sum_x, E_x);
    }

    public String var_x(HashMap<String, String> score_dict, String seq_id) {
        // Calculate and cache the variance of one sequence."
        
        if(var_cache.containsKey(seq_id))
            return var_cache.get(seq_id);

        // Fetch the mean
        (n_x, sum_x, E_x) = E_x(score_dict, seq_id);

        scores = score_dict[seq_id][1];
        num_seqs = blast_info["num_sequences"];
        logsmin = nc_info["logsmin"];

        n_others = num_seqs - n_x;

        sum_squared = numpy.sum( (scores - E_x)**2 );
        sum_squared += ((logsmin - E_x)**2) * n_others;
        var = sum_squared / (num_seqs - 1);

        // Correct for precision issues, or if no neighbors exist
        if (var < 0 || n_x==0)
            var = 0;

        var_cache[seq_id] = var;
        return var
    }

/*    public void cache_e_var(HashMap<String, String> score_dict, String query_seqs):
        self.E_cache = {}
        self.var_cache = {}

        if query_seqs is None: query_seqs = score_dict.keys()
        for seq_id in query_seqs:
            self.E_x( score_dict, seq_id)
            self.var_x( score_dict, seq_id)
        return (self.E_cache, self.var_cache)


    public void calculate_nc(score_dict, query_seqs, target_seqs, callback,
                     fd_out, seq_id_map) {
        /*Calculate NC scores for all pairs between query_seqs and
target_seqs, using the BLAST scores in score_dict.

* score_dict should be of the form:
  {'seq_id': (numpy.array([id0, id1, ...], dtype=int32),
              numpy.array([score0, score1, ...], dtype=float64)),
   ... }

  All seq_ids are expected to be of type int32.  The arrays should be
  sorted in order of ascending id.

Optional arguments:

* query_seqs and target_seqs should be lists of int32 seq_ids.  When
  unspecified, query_seqs or target_seqs will default to all keys in
  score_dict.

* callback(seq_id_0, n_0, [(seq_id_1, n_1, n_01, nc_score), ...]) will
  be called for each query sequence.  n_0 and n_1 are the number of
  BLAST hits unique to seq_id_0 and seq_id_1, respectively.  n_01 is
  the size of the common neighborhood.

* fd_out will be written with 'seq_id_0 seq_id_1 score\n'.  This must
  be an open file object.  This function will not close the file.

* seq_id_map should be a dictionary mapping from (int32) seq_id to
  some other identifier, of any type with a __str__() method.  This
  mapping will be applied before writing to fd_out.*/
/*
        if (query_seqs == null) 
            query_seqs = score_dict.keys();

        // all sequences that will be accessed should have data cached.
        // If both the queries and targets are not specified, cache all
        // sequences.
        if (query_seqs != null && target_seqs != null)
            cache_seqs = set(query_seqs).union( target_seqs);
        else
            cache_seqs = score_dict.keys();

        //print "query_seqs:", query_seqs
        //print "target_seqs:", target_seqs
        //print "cache_seqs:", cache_seqs

        if (target_seqs != null);
            target_seqs = set( target_seqs);

        // cache all sequence expectations, variances:
        (E_cache, var_cache) = cache_e_var(score_dict, query_seqs = cache_seqs);

        print >> sys.stderr, "       # Queries:", len(query_seqs)
        print >> sys.stderr, "       # Targets:", len(target_seqs) if target_seqs is not None else len(score_dict)
        print >> sys.stderr, "    NC Threshold:", self.nc_info['nc_thresh']
        print >> sys.stderr, "            Smin:", self.nc_info['smin']
        print >> sys.stderr, "Sequence DB size:", self.blast_info['num_sequences']

        t = time.time()
        for i,seq_id_0 in enumerate(query_seqs):
            if (i > 0 and i < 1000 and i % 100==0) or (i > 0 and i % 1000==0):
                print "Queries performed: %d, Elapsed time: %s" % (
                    i, datetime.timedelta(seconds = time.time() - t))

            if not seq_id_0 in score_dict:
                print >> sys.stderr, "Query '%d' not found in score_dict" % seq_id_0
                continue
            
            (n_0, nc_hitlist) = self.calculate_nc_inner(
                score_dict, E_cache, var_cache,
                self.nc_info, self.blast_info,
                seq_id_0, 
                target_seqs)

            if fd_out is not None:
                for (seq_id_1, n_1, n_01, nc_score) in nc_hitlist:
                    print >> fd_out, "%s %s %r" % (
                        seq_id_map[seq_id_0] if seq_id_map is not None else seq_id_0,
                        seq_id_map[seq_id_1] if seq_id_map is not None else seq_id_1,
                        nc_score)

            if callback is not None:
                callback( seq_id_0, n_0, nc_hitlist)

        print "Queries performed: %d, Elapsed time: %s" % (
            i, datetime.timedelta(seconds = time.time() - t))
    }*/
}
/*
class nc_base( object):


    def __init__(self):

        // try to compile and use the C module.  Otherwise, fail back to python methods
        try:
            import java.util.HashMap;

import nchelparr
            self.calculate_nc_inner = nchelparr.calculate_nc_inner
        except:
            m = """WARNING: Neighborhood Correlation C helper not found.  Using the pure
Python implementation will be significantly (~90 times) slower.  See
the README documentation for instructions on building."""
            print >> sys.stderr, m
            
import nchelparr_python
            self.calculate_nc_inner = nchelparr_python.calculate_nc_inner
*/
